Although the above (ReWrite) example has everything clearly typed, and tail recursion keeps the stack use very low, there is a problem that the Pascal version (see below) doesn't share - in the third line of primeaux,
the parameter l is used in the condition (as well as on the right hand side), which because of unique references means that the whole list of prime found so far has to be copied each time through the loop.
The following version avoids this problem by passing the whole list to primetest, which doesn't just dispose of the primes that it has tested against so far, but keeps them and returns the complete list of results as well as the result for whether the new number is prime or not. primeaux2 is 'glue' code that takes the result from primetest and adds the new number to the list or not accordingly.
Note: that this code makes use of some of the more unusual feature of ReWrite:
- primetest returns 2 results, as an arbitrary number of results can be returned;
- the list of primes is not passed as an actual list to primetest but it turned into a lot of extra arguments, and because the argument list ends with a .b, this picks up all the primes and puts them into b. This improves the code speed slightly.
nprime[n:int] -> primeaux[{},0,n,2];
(* primeaux[primes found so far, sqrt, number found, this try] *)
Time taken on LC III: 134 ticks - better than a 50-fold improvement!
It is hoped that some future version or ReWrite will be able to do this sort of code transformation automatically as an optimisation, however that is some way off.
Just for an example of how the typing helps, the above example with all the type information removed takes much longer (so putting in all those ints is worth it):
nprime[2000] -> 17389
Time taken on LC III: 556 ticks
Here is another implementation of the same algorithm, using features new to ReWrite 0.2.5 (the unrestricted use of splice on the left)
This version keeps track of four values:
* the list of primes less than or equal to the square root of the current number (this is the list to check for divisors in),
* the list of other primes found so far,
* the number of primes found so far,
* the current number being checked.
nprime[n:int] -> primeaux[{},{},n,2];
(* primeaux[(primes found so far <= sqrt[k]), (primes found so far > sqrt[k]), number found, this try] *)
(* found as many primes as we want, so finished *)
primeaux[l1:lis,l2:lis,0,k:int] -> k-1;
(* update the list of primes <= sqrt[k] if necessary *)